Coulomb explosion of uniformly charged spheroids 
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I. INTRODUCTION 

Coulomb explosion (CE) is an ubiquitous phenomenon 
in laser-matter interaction, from laser ablation and mi- 
cromachining to particle acceleration CE is the 

dominant process of ion acceleration from a cluster irra- 
diated by an intense laser pulse in the regime of so-called 
cluster vertical ionization (CVI) [3-I3l- this regime, 
the laser pulse is intense enough to remove all electrons 
from the cluster before ion motion sets in. This kind of 
ion charge state can also be generated by intense and 
short pulses of high energy photons which have become 
possible at x-ray free electron laser In both cases, 
the ion dynamics is governed by CE. 

Spherical CE has been thoroughly investigated in the 
last years due to its importance for cluster physics. In 
the case of a uniformly charged sphere, CE is self-similar 
and can be described analytically ■ The dynamics of 
CE of a non-uniformly charged sphere is more complex 
as it involves multiple flows so that a kinetic description 
is required [1, 0] . 

In contrast, ellipsoidal and spheroidal (ellipsoidal with 
a rotational symmetry) CE has been studied in the 
context of accelerator physics, where three-dimensional 
(3D) envelope equations are widely used [13, EH, or to 
model space charge effects in laser-created dense electron 
beams [Iq . 

Spheroidal clusters have also attracted a lot of atten- 
tion as they exhibit characteristic electron momentum 
distributions and due to their optical properties 

which are of great interest in, e.g., nano-optics [isl. IToj. 
Moreover, spheroidal clusters appear as a natural can- 
didate for anisotropic ion emission from clusters un- 
der intense ultrashort laser irradiation, which has re- 
cently triggered significant interest in anisotropic clus- 
ter expansion Furthermore, it has also been shown 
that in helium embedded rare gas clusters, a spheroidal 
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nanoplasma is generated by illumination with an intense 
laser pulse giving to unusual resonant heating (2l| . 

In addition, understanding spheroidal CE is crucial in 
the context of ion acceleration from a solid target irradi- 
ated by an intense, relativistic laser pulse. Recent studies 
have shown CE of thin multispecies target as a promising 
path toward high-quality ion beams [St-Ulf. Apart from 
the possibility to use CE as the principal acceleration 
mechanism, CE of the accelerated ion bunch itself has 
been shown to play a dominant role in angular as well as 
energy dispersion of ion beams generated from laser-solid 
interaction [111, [l^] . 

In this paper, we investigate CE of an initially uni- 
formly charged spheroid. In order to derive simple es- 
timates for the particle maximum energies and char- 
acteristic explosion time, we restrict ourselves to non- 
relativistic particle velocities. We then demonstrate that, 
during CE, both the spheroidal shape and uniformity 
of the charge distribution are conserved, but with time- 
dependent aspect-ratio and charge density. Therefore CE 
of a uniformly charged spheroid can be described using 
a simple, semi-analytical model for the evolution of the 
spheroid radii. This model allows us to derive the tem- 
poral evolution of the particle energy distribution and 
maximum energies (along the spheroid principal axes) 
as a function of the spheroid initial aspect-ratio, charge 
density and total charge. Our theoretical predictions 
are then compared to molecular dynamics (MD) and 3D 
particle-in-cell (PIC) simulations. These simulation tools 
are the most widely used methods to model laser-cluster 
and laser-plasma interaction. However, they are known 
to be computationally costly, so that the results obtained 
in this paper are interesting for various applications, from 
non-spherical cluster CE to non-neutralized charged par- 
ticle beam propagation through a vacuum. 

The paper is structured as follows. Section HIl presents 
our semi-analytical model. Predictions from the model 
are then compared to both MD simulations (Sec. lIIip and 
PIC simulations fSec. llVj). Finally, Sec. IVl summarizes 
our findings. 
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II. SEMI-ANALYTICAL MODEL FOR 
COULOMB EXPLOSION OF A UNIFORMLY 
CHARGED SPHEROID 



A. General considerations on uniformly charged 
spheroids 



FIG. 1. Dependence of the shape-functions Co(q) (black 
curve), (\\{a) (blue curve) and C±{ce) (red curve) on the 
aspect-ratio a. The lower row shows a prolate (cigar-shaped) 
spheroid with a = 0.1 and an oblate (disk-shaped) spheroid 
with a = 10. 



In this first Section, we lay the basis for our simple 
model of CE of a uniformly charged spheroid. To do so, 
let us first recall the electrostatic potential at a position 
X = (a;, y, z) inside a uniformly charged ellipsoid centered 
in X = and with radii Wx, Wy and Wz along the direc- 
tions X, y and z, respectively [22| : 



(x) = Trfc pc Wy Wz 



1 - 



ds 



(1) 



where k = (47reo)^^ and pc~ Zen is the charge density 
(typically Z is the mean ion charge state and n is the ion 
density), and ip{s) = {w^ + s) {Wy + s) [w1 -t- s). 

Objects with a rotational symmetry are of particular 
importance for many applications such as cluster explo- 
sion or particle acceleration. Hence, we introduce the 
radial coordinate r ~ -^/y^ + and restrict our study to 
the case of a spheroid: Wx = w\\ and Wy = Wz = 'w± so 
that V'(s) = (w^ + s) (w^ + s)^. Finally, one obtains for 
the electrostatic potential inside the spheroid: 

where we have introduced the spheroid aspect-ratio a ~ 
w±/w\\ and: 

2 Jo [a-' + s) Vl + s 

C||(a) = ^ / f 2 I \n ^ \3/2 ' 
2 Jo [a^ + s) (I + sy/^ 

C±{a)^— / ———j==. (5) 

In what follows, we refer to these functions as shape- 
functions as they depend only on the spheroid aspect- 
ratio. Equations ([3])-([5|) are here written in their com- 
pact, integral forms. These integrals can, however, be 
expressed as functions of inverse trigonometric and hy- 
perbolic functions, see e.g. Ref. The dependencies 



of these shape-functions on a are given in Fig. [TJ We 
also give their limits for the extreme prolate (cigar-shape, 
a ^ 1) and oblate (disc-shape, a 1) spheroids, as well 
as for the sphere (a = 1): 
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Equation ([2]) illustrates the well-known result that 
the electrostatic potential inside a uniformly charged 
spheroid is a quadratic function of the space coordinates. 
As for the electric field inside the uniformly charged 
spheroid, it can be easily expressed in cylindrical coordi- 
nates: 

E(a:,r) = S||(a;)i + ^i(r)f , (6) 

where x and f are the longitudinal and radial unit vec- 
tors, respectively, and: 

Ei\ (x) = -dx (l){x, r) = Airk Cy (a) x , (7) 
E^{r) = -drcj){x,r)=A^kp,(:i_{a)r. (8) 

Interestingly, beside the non-trivial dependency on the 
aspect-ratio a, the longitudinal componant of the elec- 
tric field inside the spheroid is a linear function of x only, 
while the transverse componant is a function of r only. 
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The non-relativistic equations of motion in such an elec- 
tric field for a particle with charge Ze, mass m and initial 
position (xo,ro) simply read: 



dt^ 



—x = KCi\{a)x, 
f = KCj_(a)r, 



(9) 
(10) 



where a and k ~ Ank pc (Ze/m) depend only on time, 
and X ~ x/xq and r = r/rg. Considering initial con- 
ditions x\t=a = f\t=o = 1 and dtx\t=o = dtf\t=o = 0, 
Eqs. © and pU)) are found to be independent on the ini- 
tial coordinates xq and tq. As a consequence, a particle 
initially located at a position {xo,ro) will subsequently be 
at a position (xq x,rof) where x and f do not dependent 
on the initial position. Therefore, it is straightforward 
to obtain that (i) CE of an initially uniformly charged 
spheroid conserves the spheroidal shape (albeit, as we 
will see in Sec. Ill B[ with a time-dependent aspect-ratio), 
and that (ii) the charge distribution inside the spheroid 
remains uniform (albeit time-dependent). 

Although the above calculations have been performed 
considering Coulomb interaction in a spheroid, we would 
like to stress that similar conclusions can be drawn for 
the more general case of an ellipsoid, as well as for any 
Coulomb-like force, such as, e.g., gravitation [23 |. 

The above considerations strongly simplify the mod- 
elling of CE of a uniformly charged spheroid. The prob- 
lem can now be solved by considering equations for the 
evolution of the longitudinal and transverse radii of the 
spheroid. 



B. Coulomb explosion of a uniformly charged 
spheroid 

1. Governing equations 

Let us consider a uniformly charged spheroid with ini- 
tial radii w^^^q and w_l,o, initial ion charge density uq. 
Obviously, the total charge Q = (47r/3) ifn^o ai^ e tiq) 
is conserved during explosion. The non-relativistic equa- 
tions of motion ^ and pi7)) are also valid for particles 
initially located on the outer shell of the spheroid at 
{x = w\\fl,r = 0) and (x = 0,r = w^^q). Then, using 
Eqs. ([7]) and it is straightforward to derive a system 
of two second order differential equations on the time- 
dependent longitudinal and transverse radii w\\ and w^: 



df' 

dt^' 



WW 



^"±.0^11,0 



w± 

Wll 



2 ^±.0«^ll,0 . fwi_ 

^po '■ C-1 — 



where we have introduced the plasma frequency: 



ujpo = ^/Z^e^ no/(eo™) 



(11) 
(12) 

(13) 



Let us now normalize the time to uj^^^ (r = ojpo t) , the 
longitudinal radius to (^n ~ w\\/w\\,i^) and the trans- 
verse one to wj^.o {w± = w^/w±^q). The system of 
Eqs. (HH) and ^ then reads: 



dr2 



— C|| "0 — 

^1 V ^11 



1 



w± 
"0 — 

W\\ 



(14) 
(15) 



where ao = w^J-.o/wy^o is the spheroid initial aspect-ratio. 
Assuming that all particles in the spheroid have initially 
no velocity, one has for initial conditions: 



dr 



W||(r = 0) = 1,u;_l(t = 0) = 1, (16) 
0)=0. (17) 



d , , d 
— W||(r = 0) = , —w I (T 



dr 



Note that Eq. (IT71) also implies that particles have no ini- 
tial temperature. For bare ion spheroids, this situation 
arises when electron removal is fast enough for ion heat- 
ing through electron-ion collisions to be negligible. Under 
such circonstances, the hypothesis of uniform charge den- 
sity should also be verified as long as the initial atomic 
density is uniform. 

Before discussing in more details the solution of this 
system, we note that, using these normalizations, veloci- 
ties in the longitudinal and transverse directions are nat- 
urally expressed in units of ^^^11,0 and w±,07 respec- 
tively. Correspondingly, energies in the longitudinal and 
transverse directions are normalized to q ~ Sq/olq and 
^±,0 = OLQ^o, respectively, where we have introduced the 
characteristic energy: 



(18) 



Due to the complex dependency of the shape- functions 
on the time-dependent aspect-ratio, no general (for any 
ao) analytical solution can be obtained and the sys- 
tem of Eqs. (HH) and ([T5|) has to be solved numerically. 
However, analytical solutions of the system can be ob- 
tained in the case of the sphere (where the aspect-ratio 
remains constant in time w^l/ww ~ ~ 1), in the 
case of an infinitely large disc (ao — > +oo and consid- 
ering uqw^/ww — > +00 for all times), and in the case 
of an infinitely long cylinder (ao — >■ and considering 
cko w^_/'W\\ — ?> for all times). We present briefly the an- 
alytical solutions for these particular cases (Sec. IIIB 2[) 
before discussing in more details the numerical solutions 
for arbitrary initial values of ao (Sec. IIIB"3| . 



2. Particular cases 

a. Coulomb explosion of a uniformly charged sphere 
Spherical CE (ao = 1) has been widely studied in the 
context of many applications, e.g. in cluster physics 0- 
01 • In this case, the system of Eqs. and reduces 
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to a single differential equation on the sphere radius R 
WW = w±: 



dr' 



1 



3i?2 



(19) 



The first integral of Eq. (jl9p is obtained after multiplying 
both sides by -^R and integrating from to r. Consider- 
ing that all particles have initially no velocity, we obtain: 



(20) 



This equation describes the transformation of potential 
energy [right-hand-side (rhs) of Eq. ([20| ] to kinetic en- 
ergy [left-hand-side (Ihs) of Eq. ((20|) ] for a particle lo- 
cated on the outer shell of the spheroid. In our normal- 
ized units, energies £ arc expressed in units of g = 
£±fl ^ £q = mujpQRl^, where Rq is the initial radius of 
the sphere and ujpQ is given by Eq. (|13p . The final kinetic 
energy of an ion of the sphere outer shell is therefore: 



(21) 



The autonomous differential equation (|20p has a formal 
implicit solution: 



dr 



r- 1 



(22) 



This solution and the temporal evolution of the outer 
shell kinetic energy Smax a-re shown in Fig. [21 On long 
time scales (r = ujpot > 10), we find R ^ \f2j2>T^ i.e., 
most of the potential energy has been transformed into 
kinetic energy and the sphere expands with the constant 
velocity -^2/3 w^o -Ro- On shorter time scales we actually 
observe the Coulomb explosion, i.e., 50% of the poten- 
tial energy is transformed in kinetic energy after a time 
r ~ 2.8, while 80% is transformed after r ~ 7.2. The 
characteristic time scale of spherical explosion is there- 
fore the inverse initial plasma frequency . 

Following Refs. 0, H, we can derive an analytical 
expression for the asymptotic (t — > oo) particle energy 
distribution. As previously underlined, in our model, the 
electric field inside the sphere is a linear function of the 
radius and particles do not overtake each other during 
expansion. The electric field seen by a particle initially 
located at f g < 1 can thus be easily obtained as a func- 
tion of the charge (;(ro) inside the sphere with normalized 
radius fo. The equation of motion for this particle then 
reads: 



'0 

3f2 



(23) 



Then the first integral simply reads (assuming zero initial 
velocity): 



dT 



(24) 



This equation once more describes the energy conserva- 
tion and shows that an ion initially located at a position 
fo has obtained, at the end of the acceleration process, a 
kinetic energy f (fp) = r^jZ. Now, the normalized radial 
particle density at initial position fo is simply: 



dN 
dro 



3f,2 



OH 



(25) 



where 9h is the Heaviside function, from which we derive 
the ion energy distribution at i — > oo: 



dN _ dN dfo 
dS dfo d£ 



^3 5 6*^(1/3 -£) 



(26) 



Hence, one obtains that the particle energy distribution 
scales as the square-root of the ion energy up to the max- 
imum energy £°° = 1/3. 

At this point, we want to stress that this asymp- 
totic energy distribution can actually be generalized to 
all times t. It is indeed well-known that CE of a uni- 
formly charged sphere is self-similar and that the veloc- 
ity distribution inside the sphere increases linearly with 
the distance to its center. From this we can derive the 
fraction of particles with an energy below £ < £max{t)- 
N{£) — [£ /£max{t)]'^^^ , where £max{t) is the maximum 
particle energy at time t. We finally obtain the time- 
dependent spectrum: 



dN 



V£ 



2 £-^/2 (t) 



^H[^rnax{t) — S] 



(27) 



h. Coulomb explosion of a uniformly charged, in- 
finitely long cylinder In the case of a uniformly charged, 
infinitely long cylinder, the system of Eqs. p4)) and 
reduces to a single differential equation for the cylinder 
radius wj,: 



1 



dT^^ 



2iD_L 



(28) 



Once more, the first integral of Eq. ((28| describes energy 
conservation: 



i In 



(29) 



In contrast to the spherical case considered above, the 
logarithmic potential on the rhs of Eq. goes to infin- 
ity for increasing w^. This unphysical behavior follows 
from our choice of an infinitely long (thus with infinite to- 
tal charge) cylinder. As a result, the energy of the outer 
shell artificially formally diverges. 

Again a formal implicit solution can be obtained for 
the cylinder radius: 



dw 



(30) 



The temporal evolution of zDj^ and the outer shell energy 
f _L are presented in Fig. [2] Expansion occurs with an in- 
creasing velocity, and no saturation of the kinetic energy 
is observed. 
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FIG. 2. Temporal evolution of a) the normalized spheroid 
radii {R/Ro, W||/u;||_o and w±/'W±fi), and b) the maximum 
energies {£max/£o, £\\,max/£\\,o and £±,max/£±,o) for spherical 
explosion (solid curves), cylindrical explosion (dashed curves) 
and planar explosion (dot-dashed curves), respectively. 



c. Coulomb explosion of a uniformly charged, in- 
finitely large disc In the case of a uniformly charged, 
infinitely large disc, the system of Eqs. (fT4|) and (fT5|) re- 
duces to a single differential equation for the disc thick- 
ness : 

|^^„=1. (31) 

In this case, the disc thickness increases due to a constant 
electrostatic field. It simply reads: 

*|| = Y + 1 • (32) 

This result is shown in Fig. [2] In this case as well, no sta- 
tionary state is obtained and the kinetic energy increases 
arbitrarily. 

3. Numerical solutions for an arbitrary initial aspect-ratio 

As previously discussed, in the general case (for any 
initial aspect-ratio ao), the system of Eqs. (HH) and 
must be solved numerically. In this paper, it is done using 
a simple Eulcr method. Numerical solutions for different 
values of ao are now discussed. 

Figure |3] shows the temporal evolution of the longitudi- 
nal and transverse radii of the spheroid (Figs. [5^ and[5]3, 
respectively) and of the longitudinal and transverse ki- 
netic energies (Figs. [31; and[3ji, respectively). It is com- 
plemented by Fig. 2] which shows, as a function of the ini- 
tial aspect-ratio ao, the times required for the kinetic en- 
ergy (for purely longitudinal or purely transverse motion) 




FIG. 3. Temporal evolution of spheroidal Coulomb explosion: 
a) longitudinal radius of the spheroid; b) transverse radius of 
the spheroid; c) energy associated with longitudinal motion; 
and d) energy associated with transverse motion. Color codes 
are as follows: Prolate spheroid (dashed curves, ao ~ W^^ 
red, ceo ~ green and qq = 10~^ blue). Sphere (solid 

black curve, qo = 1). Oblate spheroid (solid curves, ao = 10'^ 
red, Qfo = 10^ green and qq ~ W blue). The grey solid curves 
in panels a) and b) show analytical predictions for planar 
(ao — >■ oo) and cylindrical (a — >■ 0) explosion, respectively. 
Note that energies in panels c) and d) are normalized to their 
asymptotic values at t — >■ cx3. 



to reach 50 % or 80 % of its maximum value (Fig. 
the corresponding aspect-ratio of the spheroid at these 
times (Fig. |4}d), the final aspect-ratio aoo (Fig. Ht) and 
the final energies normalized to £|| q and £±,o (Fig. HJi). 

Let us first address the case of spherical explosion 
(ao = !)• Numerical solutions allow us to recover the 
analytical findings of Sec. Ill B 21 The radius and energy 
evolutions (Figs. [3]) are in perfect agreement with what is 
presented in Fig.[5J It is also worth pointing out that for a 
given initial charge density (i.e., for a fixed value of Wpo), 
spherical expansion is faster (in terms of energy conver- 
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FIG. 4. Dependence on the initial aspect-ratio qq of: a) the 
time to reach 50 % (dashed curves) or 80 % (solid curves) of 
the maximum kinetic energy; b) the corresponding intermedi- 
ate aspect-ratio qq iSx/wn ; c) the final aspect-ratio aoo', d) the 
final energies £^ and £]° normalized to £\\fi and £±.o, respec- 
tively. In panels a), b) and d), blue curves account for motion 
in the longitudinal direction, and red curves for motion in the 
transverse direction. 



sion) than for any other initial aspect-ratio oq. This can 
be observed in Figs.[3t and Figs.[3ji, as well as in Fig.|4^, 
where the characteristic time to reach a given fraction of 
the final kinetic energy reaches a minimum for ag = 1. 
Finally, we recall that the characteristic time for spheri- 
cal CE is of the order of uj~q and the final velocities of the 
outer shell in the longitudinal and transverse directions 
are W||,max(i oo) = wj_,max(i ^^00) = ^/2J3uJpO Rq. 

In the case of a prolate (cigar-shape, ao ^ 1) spheroid, 
expansion occurs mainly in the transverse plane. This in- 
tuitive result is illustrated in Figs. [3^ and [3)3 where the 
transverse radius of the spheroid increases much faster 
than the longitudinal one. Note that, even though the 



transverse radius of the prolate spheroid evolves initially 
faster than in the spherical case, spherical expansion re- 
mains faster in terms of conversion from potential to ki- 
netic energy (Fig. |3Ji) . The quasi-stationary state of ex- 
pansion (i.e., expansion at quasi-constant velocity) is in- 
deed reached later for ao 7^ 1. Note also that saturation 
in the kinetic energy arises once the spheroid assumes a 
quasi-spherical shape: when the kinetic energy reaches 
80 % of its maximum value, the spheroid aspect-ratio is 
indeed quite close to unity (Fig. |3}d) . The final aspect- 
ratio is nevertheless much larger than unity as the final 
transverse velocity is much larger than the longitudinal 
one. We also see from Fig. that spheroidal expansion 
in the limit ao ^ 1 occurs on a time scale larger than 
(aQLUpo)^^ . Similarly, we see in Fig. HJi that the asymp- 
totic {t ^ 00) energy for purely longitudinal motion £^ 
does not exceed a^fn q = ao£o while the corresponding 
energy for purely transverse motion £^ > fo/S (for ex- 
ample, for 10~^ < ao < 10~^, we find that £]° ranges 
between £±,0 = O!o£o and 7£±,o ~ 7ao£o)- The final 
transverse energy is therefore significantly larger than 
the longitudinal one (see also theoretical predictions in 
Fig. [5]). This is a consequence of the initial geometry, and 
it is responsible for the final oblate shape of the spheroid 
observed in Fig. HJd. 

Let us now focus on the case of an oblate (disc-shape, 
ao ^ 1) spheroid which is particularly interesting when 
considering laser-generated ion bunches from a solid tar- 
get. As expected, expansion initially occurs in the lon- 
gitudinal direction (Figs. andOa). Transverse expan- 
sion eventually occurs later, once the spheroid longitu- 
dinal radius becomes comparable to its transverse one. 
At this time, the spheroid aspect-ratio becomes close to 
unity and a non-negligible fraction of the potential energy 
has already been converted into kinetic energy (Fig.Hb). 
Finally, one can extract from the numerical results the 
characteristic expansion time in the limit of large ini- 
tial aspect-ratio to be oc ^Jao/ujpo- The asymptotic 
maximum energies can also be easily extracted: £^ ~ 

0.63ao'S'||,o = 0.63 £0 and £^ ~ 0.46£_L,o/ao = 0.46 fo- 
We thus obtain that the final energies in the longitudi- 
nal and transverse directions arc of the same order. This 
leads to a final aspect-ratio aoo ^ ^J^i^ ~ ^'^^ close 
to unity, as shown in Fig. H};. 



C. Energy spectra 

Our model for spheroidal CE allows us to derive the 
maximum energies for motion along the spheroid prin- 
cipal axes at time t. As the density inside the spheroid 
remains uniform, the velocity distribution along the prin- 
cipal axes has to be a linear function of spatial coordi- 
nates: 

V\\{x) = {x/w^)v\\^rnax , (33) 

v±{r) = {r/un_) v^.max , (34) 
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FIG. 5. Schematic representation of equivelocity surfaces at 
fixed time t during CE of: (left) a prolate (cigar-shaped, 
ao < 1) spheroid, and (right) an oblate (disk-shaped, ao > 1) 
spheroid. In this two-dimensional representation, equiveloc- 
ity surfaces correspond to concentric ellipses with the same 
aspect-ratio. 



where w± are the spheroid radii, and v\\^rnax ^iid 
v±.max the particle maximum velocities at a given 
time t. Considering an homogeneous charge density in 



III. MOLECULAR DYNAMIC SIMULATIONS 



Molecular dynamics [2^ simulations of CE of an (ini- 
tially) uniformly charged spheroid arc now discussed. To 
initialize our simulations, N = 5000 particles [here, hy- 
drogen ions {Z = 1, m = 1836 mg, where me is the 
electron mass)] are randomly placed within a spheroidal 
volume so that the initial particle density inside this vol- 
ume is homogeneous. Here we chose an atomic density 
no ~ 9.7 X 10^^ cm""^ (correspondingly, the sphere ra- 
dius is Rq ~ 2.3 nm), which is characteristic of hydro- 
gen clusters. To avoid unphysically large contributions 
to the energy spectrum, we enforce a minimum inter- 
particle distance (~ 75 % of the average interparticle 

— 1/3 

distance dmm ^ n-Q ). Furthermore, all particles were 
taken initially at rest. Then, knowing the initial state of 
all particles, we solve Newton's (non-relativistic) equa- 
tions of motion for each of them using the velocity Vcrlct 
scheme [13, [13 and direct calculation of the Coulomb 



the spheroid, one can easily derive the time-dependent 
normalized energy distribution for motion along the lon- 
gitudinal and transverse directions [cf. Eq. (P?]) ]: 

dN 3 7^ 

^=2^3/2 ^H{t\\.,nax-t-\\), (35) 
II ,7nax 

dN 3 vn 

where £\y^ax = "^^^maxl^ aild ■?J_,max = mv\,^^j2 

are time-dependent and derived from our model. 

The total energy spectrum can also be derived by 
considering equivelocity surfaces as concentric homeoids 
(spheroidal surfaces). We obtain from Eqs. ([33| and (p4|) 
that these homeoids are actually similar, i.e. they have 
the same aspect-ratio ct,, = {w\_lw\^v\^^ma,xl'v^,max at 
fixed time t (see Fig. [5]) [l^] . This allows us to calculate, 
for a given energy f , the fraction A^(f ) of particles in 
the spheroid with a lower energy, and finally derive the 
normalized energy distribution. 



(37) 



(38) 



forces between all ions. 

Several simulations were performed only changing the 
initial aspect-ratio in the range ao = 0.1 — 10. Figure [S] 
shows the maximum energy for longitudinal and trans- 
verse motion as predicted by our semi-analytical model 
(Sec. ini) and as extracted from MD simulations. We 
stress that, here, energies are normalized to the maxi- 
mum energy £5 [Eq. (PT|)] resulting from CE of a sphere 
with similar density and total charge (in practical units, 
£5 ~ 3.1 kcV under current conditions). Figure [H] shows 
a rather good agreement between our simplified model 
[solutions of Eqs. and (|15p ] and simulations. Also 
note that MD results confirm the theoretical prediction 
(clearly shown in Fig.[S]) that, for a given total charge and 
charge density in the spheroid, the maximum longitudi- 
nal (transverse) energy is obtained for a slightly oblate 
(prolate) spheroid. 

To further investigate MD simulation results, we 
present in Fig. [7] the energy spectra obtained at the 
end of the simulation [at a time t ~ 400 ijJ~^ (and 



For a prolate (cigar-shaped, ao < 1) spheroid, the total energy spectrum reads: 



d£ 2 I 1 ^ / £j_,„iax-£ r c < f < f , 

For an oblate (disk-shaped, ao > 1) spheroid, the total energy spectrum reads: 



In the next Sees. Illll and llVj we compare these predictions from our model to MD and PIC simulations. 

I 



8 




FIG. 6. Dependence of the maximum energy for the longi- 
tudinal and transverse motions (blue and red curves, respec- 
tively) on the initial aspect-ratio qq. Solid lines correspond to 
semi-analytical predictions [Eqs. (|14|) and (|15|) ]. circles to MD 
simulation results, and triangles to PIC simulation results ob- 
tained at the end of the simulation. Note that energies are 
normalized to the asymptotic energy £s = fo/3 [Eq. (|21|) ] 
obtained for spherical explosion. 
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FIG. 7. Energy spectra obtained from MD simulations at 
t ~ 400 i^j7o^: a,b) for ao = 1, c,d) for qq ~ 0.1, and e,f) for 
ao = 10. The left panels a,c,e) show the directional spec- 
tra: the blue curves correspond to particles emitted within 
an angle ±7r/20 of the longitudinal {x) direction, the green 
dashed curves and red curves to particles emitted within an 
angle ±7r/20 of the transverse {y and z) directions, respec- 
tively. The right panels b,d,f) show the total spectra (all 
particles are accounted for). The gray dashed lines show the 
theoretical predictions from our model. 



cjpQ^ ~ 2.4 fs)] and compare them to theoretical predic- 
tions from Eqs. (P?]) and ([55]) -([55 ]) . Distributions in lon- 
gitudinal (transverse) energy are obtained by consider- 
ing particles emitted within an angle 9 = ±7r/20 around 
the longitudinal (transverse) direction. Panels a) and b) 
correspond to spherical CE (ao = 1). In this case, the 



energy distribution is in very good agreement with the 
theoretical prediction from Eq. (P7)) for energies up to 
90 % of the maximum energy. Near the maximum en- 
ergy, however, a peak appears in the simulation which is 
not predicted by our model. 

Figures [Thjd and[7^,f show similar results for CE of a 
prolate (cigar-shape, ao = 0.1) spheroid and an oblate 
(disk-shape, ao = 10) spheroid, respectively. Once more, 
a very good agreement is obtained between theoretical 
predictions from our model and MD simulations. Only 
at the maximum energy a peak is observed for any direc- 
tion of emission in the MD spectra, which is absent in the 
model spectra. This peak is a consequence of the discrete 
nature of particles, which leads to a gradual decrease in 
the particle density at the surface of the spheroid. The 
thickness of the corresponding surface layer is of the or- 
der of the average intcrparticle separation. The decreas- 
ing density results in ion wave breaking (or formation of 
a shock), whose characteristic signature is a peak in the 
energy spectrum BSHi,!!!. The peak is found to be 
sensitive to the initial particle distribution and its mag- 
nitude decreases with the increase of the total number 
of particles. This behavior is captured by the MD cal- 
culations, which take into account the motion of discrete 
ions, but is neglected in our simplified model, which con- 
siders the evolution of a continuous particle density with 
a sharp cutoff at the surface. 



IV. PARTICLE-IN-CELL SIMULATIONS 

We now present results from simulations of spheroidal 
CE obtained using the massively parallel 3D PIC code 
CALDER [3l|. The PIC simulation technique consists 
in solving the Maxwell- Vlasov system, and thus offers 
a mean-field kinetic description for the plasma dynam- 
ics [s^]. In PIC codes, the Vlasov equation is solved 
by discretizing the particle distribution functions as a 
sum of so-called macro-particles and by solving, for 
each of these macro-particles, the corresponding (rela- 
tivistic) equation of motion in the electromagnetic field. 
Then, the Maxwell- Ampere and Maxwell-Faraday equa- 
tions are solved on a Yee-mesh using the finite-difference 
time-domain method (ssj . This numerical scheme, cou- 
pled to the standard current and charge deposition algo- 
rithm in a PIC code, does not automatically satisfy the 
Poisson equation, which has to be enforced by correct- 
ing the electric fields at each time step. In CALDER, 
this is done by using the usual technique proposed by 
Boris [13]. This study presents the first simulations per- 
formed with CALDER in the case of an initially strongly 
non-neutral plasma for which it is of the utmost impor- 
tance to accurately correct the electric field at all time- 
steps. This difficulty can however be alleviated by using 
charge-conser ving algorithms such as the one proposed 
by Esirkcpov [35[ . In this case indeed, Poisson equation 
has to be solved only at the first time-step. 

Now, we present simulation results of spheroidal CE for 
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FIG. 8. Results from PIC simulations for ao = 1. a) Tem- 
poral evolution of the sphere radii, b) Temporal evolution 
of the maximum kinetic energies, c) Ion energy spectrum at 
t ~ 10.6 i^^Q ■ Quantities in panels a) and b) are presented 
for all three spatial directions: along the x-direction (blue), 
y-direction (dashed green) and z-direction (red). The gray 
dashed lines show theoretical predictions from our model. 
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three different values of the initial aspect-ratio ag = 0.1, 
ao = 1 and ao = 10. All three spheroids consist of fully 
ionized carbon ions (Z = 6, m = 12 x 1836 me) with 
density no ~ 9.2 x 10^^ cm~'^ and total charge (5—19 pC. 

We first consider CE of a sphere (ao = 1) with ini- 
tial radius i?o — 80 nm. In this simulation, the mesh 
size in all three directions is as small as Ace = Ay = 
Az = i?o/20, and 200 macro-particles per cell are used. 
Figures and [SJ) display the temporal evolution of the 
spheroid radii (normalized to their initial value i?o) and 
of the maximum kinetic energy [normalized to the fi- 
nal energy predicted from the model Es (in practical 
units Es — 12-8 MeV for present parameters)], respec- 
tively. Note that these quantities have been extracted 
along all three spatial directions (a;,y and z), confirming 
that the CE dynamics remains spherical within a ~ 2 % 
error. Furthermore, an excellent agreement is found be- 
tween our theoretical model and simulations concerning 
the temporal evolution of both the maximum energy and 
sphere radius. Figure [5J; shows the carbon energy dis- 
tribution at time t ~ 10.6 wIq' (for current parameters 



FIG. 9. Results from PIC simulations for ao = 0.1. Tem- 
poral evolution of the spheroid radii: a) in the longitudinal 
direction, b) in the transverse direction, c) Temporal evolu- 
tion of the maximum energies, d) Corresponding spectra at 
t ~ 45.5 ij-ipo^, the total spectrum is also shown as a black solid 
line. Colors codes are chosen as in Fig. [T] The gray dashed 
lines show theoretical predictions from our model. 



u)^Q ~ 4.6 fs) and their comparison to theoretical pre- 
diction from Eq. (j27|) . Here again a very good agree- 
ment is found between our theoretical model and PIC 
simulations. It is interesting to see that, in these sim- 
ulations also, a peak is present at maximum energy. It 
also originates from a smoothly decreasing ion density at 
the sphere edge leading to the formation of a shock. In 
contrast to MD simulations, however, here the decreasing 
density region is not the result of the discrete particle dis- 
tribution, but of the projection of the particle density on 
the mesh. Note also that, at the end of the simulation (at 
t ^ 10.6 (^pq), the maximum ion energy reaches ~ 85 % 
of its theoretical maximum value (for cUpQ t — >■ oo), which 
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FIG. 10. Results from PIC simulations for ao — 10. Tem- 
poral evolution of the spheroid radii: a) in the longitudinal 
direction, b) in the transverse direction, c) Temporal evolu- 
tion of the maximum energies, d) Corresponding spectra at 
t ~ 33.3 , the total spectrum is also shown as a black solid 
line. Colors codes are chosen as in Fig. [T] The gray dashed 
lines show theoretical predictions from our model. 



is in good agreement with the dynamics illustrated in 
Figs.[l;,d. 

Let us now consider the case of a prolate (cigar-shape, 
aQ <^ 1) spheroid with longitudinal radius w\\ q = 4.6 Rq 
and transverse radii wj_^o ~ 0.46 Rq (corresponding to 
an initial aspect-ratio ao = 0.1). In this simulation, the 
mesh size is Aa; = Ay = Az = i?o/10 and 500 particles 
per ceh have been used. Figures inti,b,c display the tem- 
poral evolution of the normalized spheroid radii and max- 
imum kinetic energies along all three space dimensions, 
respectively. Figure IHli shows the energy distributions at 
the end of the simulation together with their comparison 
to theoretical predictions from Eqs. (|35p . (pS)) and ([571) . 



A fair agreement is found between PIC simulations and 
analytical predictions. 

Finally, the case of an oblate (disk-shape, ap ^ 1) 
spheroid with ao = 10 (w||,o = ^o/5 and w±^o = 2i?o) is 
presented in Fig. 1101 The mesh sizes in this simulation 
were set to Ax = Ra/20 and Ay = Az = i?o/10 and 
each cell initially contained 300 particles. The temporal 
evolution of the normalized radii and normalized maxi- 
mum kinetic energy along the three spatial directions are 
displayed in Fig. [TOb.b andfTOb. respectively. Figure [TOH 
shows the energy distributions and their comparison to 
theoretical predictions from Eqs. ([55]) . ([55]) and ([55)) . A 
fair agreement between the PIC simulations and the an- 
alytical results is also obtained in this case. 

These simulations demonstrate that our simple model 
correctly describes the CE dynamics of an initially uni- 
formly charged spheroid. We attribute discrepancies be- 
tween PIC simulation and our model predictions to the 
limited resolution of the numerical mesh. Due to tech- 
nical constraints of our computing facilities, we are cur- 
rently not able to run simulations with a higher resolu- 
tion. 



V. CONCLUSION 

We have developed a simple, semi-analytical model for 
CE of a uniformly charged spheroid. In the limit of non- 
relativistic particle velocities, this model gives access to 
the maximum energy a particle can reach at a given time, 
the time-dependent particle energy distributions, and the 
characteristic time of CE. All these quantities can be de- 
fined as a function of the spheroid aspect-ratio, charge 
density and total charge. Our theoretical predictions are 
found to be in remarkably good agreement with particle 
(both MD and PIC) simulations. 

As 3D kinetic simulations come at a high computa- 
tional cost, our results are particularly useful when con- 
sidering acceleration of ions in the pure CE regime, origi- 
nating from (spherical or non-spherical) clusters, or from 
thin, solid targets. 

Indeed, our results should be directly applicable in the 
so-called CVI regime. This regime where electrons are 
expelled from the cluster on a time much shorter than 
the characteristic time of ion motion can be accessed by 
using either ultra-intense laser or x-ray pulses. 
Moreover, with the recent progress in nanotechnology, 
CE of nanostructured targets can be considered. Our re- 
sults may thus give us simple design guidelines how to 
optimize target properties, e.g., for inertial fusion appli- 
cations [3^ or to maximize ion collision events for neu- 
tron production [33l . 

Last but not least, our results can also be helpful to 
model laser-solid target interaction for ion acceleration 
which is characterized by the emission of short, compact, 
and highly charged ion bunches. Propagation of these 
bunches, e.g. throu gh a vacuum, is strongly affected by 
space charge effects fill [l^ . 



By approximating the accel- 
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erated ion bunches as uniformly charged spheroids, the CCRT and CINES (Grant 2010-x2010056304). 
results presented here may allow us to derive the condi- 
tions required for limited energy and angular dispersions. 
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